 data {
     real low;
     real high;
     int<lower=0> N;
     int<lower=0> K1;
     int<lower=0> K2;
     int<lower=1> C; #n counties
     int<lower=1> P; #n parties
     int<lower=1, upper=C> county[N];
     int<lower=1, upper=P> party[N];
     int<lower=1>N_val;
     vector<lower=low, upper=high>[N_val] val;
     matrix[N,K1] x1;
     matrix[N,K2] x2;
     int<lower=0> y[N];
 }

 parameters {
     vector[K1] beta1;
     vector[K2] beta2;
     vector[C] ranef_c_1;
     vector[C] ranef_c_2;
     vector[P] ranef_p_1;
     vector[P] ranef_p_2;
     real alpha1;
     real alpha2;
     real<lower=0, upper=100> sigma_b1;
     real<lower=0, upper=100> sigma_b2;
     real<lower=0, upper=100> sigma_p1;
     real<lower=0, upper=100> sigma_p2;
     real<lower=0, upper=100> sigma_c1;
     real<lower=0, upper=100> sigma_c2;
     real<lower=3, upper=7> nu;


 }

 transformed parameters {
     vector[N] eta1;
     vector[N] eta2;
     vector[N] y_hat1;
     vector[N] y_hat2;
     eta1 = x1*beta1;
     eta2 = x2*beta2;

     for(n in 1:N){
        y_hat1[n]=alpha1+eta1[n]+ranef_p_1[party[n]]+ranef_c_1[county[n]];
        y_hat2[n]=alpha2+eta2[n]+ranef_p_2[party[n]]+ranef_c_2[county[n]];
     }

 }

 model {
    nu~gamma(4, 0.1);
    sigma_b1~normal(0,1);
    sigma_b2~normal(0,1);
    sigma_p1~inv_gamma(5,5);
    sigma_p2~inv_gamma(5,5);
    sigma_c1~inv_gamma(5,5);
    sigma_c2~inv_gamma(5,5);
    alpha1~normal(0,0.5);
    alpha2~normal(0,0.5);
    to_vector(beta1)~student_t(nu,0, sigma_b1);
    to_vector(beta2)~normal(0, sigma_b2);
    to_vector(ranef_c_1)~normal(0,sigma_c1);
    to_vector(ranef_c_2)~normal(0,sigma_c2);
    to_vector(ranef_p_1)~normal(0,sigma_p1);
    to_vector(ranef_p_2)~normal(0,sigma_p2);

     for (i in 1:N) {
         if (y[i] == 0) {
             /* Zero case */
             target += log_sum_exp(/* Structural zero */
                                   bernoulli_logit_lpmf(1 | eta1[i]),
                                   /* Poisson zero */
                                   bernoulli_logit_lpmf(0 | eta1[i]) +
                                   poisson_log_lpmf(0 | eta2[i]));
         } else {
             /* Non-zero case */
             /* First term means not structural zero. */
             target += bernoulli_logit_lpmf(0 | eta1[i]) +
                 /* y[i] is relevant only here. */
                 poisson_log_lpmf(y[i] | eta2[i]);
         }
     }
 }

 generated quantities {
     int y_rep[N];
     vector[N] log_lik;
     vector[K2] fd_beta2;
     vector[K1] fd_beta1;
     vector[N_val] pred_val;


     for (i in 1:N) {

         if (bernoulli_logit_rng(eta1[i]) == 1) {
             /* Structural zero */
             y_rep[i] = 0;
         } else {
             /* Not structural zero */
             if (eta2[i] > 20) {
                 /* To avoid erros like the below during the warmup. */
                 /* Check posterior predictive. */
                 y_rep[i] = poisson_log_rng(20);
             } else {
                 y_rep[i] = poisson_log_rng(eta2[i]);
             }
         }


         if (y[i] == 0) {
             /* Zero case */
             log_lik[i] = log_sum_exp(/* Structural zero */
                                      bernoulli_logit_lpmf(1 | eta1[i]),
                                      /* Poisson zero */
                                      bernoulli_logit_lpmf(0 | eta1[i]) +
                                      poisson_log_lpmf(0 | eta2[i]));
         } else {
             /* Non-zero case */
             /* First term means not structural zero. */
             log_lik[i] = bernoulli_logit_lpmf(0 | eta1[i]) +
                 /* y[i] is relevant only here. */
                 poisson_log_lpmf(y[i] | eta2[i]);
         }

     }
     // nationla party, beta5 mode is 1; left_party, beta7, mode is 1
     for(k1 in 1:K1){
       if (k1<5)
         fd_beta1[k1]=(exp(alpha1+beta1[5]+beta1[7]+beta1[k1]*(mean(x1[,k1]+sd(x1[,k1]))))/(1+exp(alpha1+beta1[5]+beta1[7]+beta1[k1]*(mean(x1[,k1]+sd(x1[,k1]))))))-(exp(alpha1+beta1[4]+beta1[6]+beta1[k1]*(mean(x1[,k1]+sd(x1[,k1]))))/(1+exp(alpha1+beta1[4]+beta1[6]+beta1[k1]*(mean(x1[,k1]+sd(x1[,k1]))))));
       else
           fd_beta1[k1]=exp(beta1[k1])/exp(1+(beta1[k1]));
 }
      for(k2 in 1:K2){
         if (k2<5)
          fd_beta2[k2]=exp(alpha2+beta2[5]+beta2[7] + beta2[k2]*(mean(x2[,k2]+sd(x2[,k2]))))-exp(alpha2+beta2[4]+beta2[7] + beta2[k2]*(mean(x2[,k2]-sd(x2[,k2]))));
        else
           fd_beta2[k2]=exp(alpha2+beta2[k2]*1)-exp(alpha2+beta2[k2]*0);

      }

      for(n in 1:N_val){
        pred_val[n]=exp(alpha2+beta2[5]+beta2[7]+beta2[1]*val[n]);
      }


 }
